##=============================================================================#
## Code to replicate the figures for the "Opposition Parties and the Timing of No-Confidence Motions", PSRM
## Based on Daina Chiba's code
## Created : 1-21-11
## Modified: 12-21-14
##=============================================================================#

rm(list=ls())
require(lattice)
require(MASS)
require(ggplot2)
require(gridBase)
require(gridExtra)
library(foreign)

# Set up working directory
#setwd("")

###############################################################################################
## Figure 2: Changes in the predicted probability of opposition proposing a no-confidence motion over the baseline probability
###############################################################################################
b <- read.csv("prob_p.csv", sep=",", header=TRUE)

mat <- matrix(NA,9,ncol=6)         
colnames(mat) <- c("pr","lower.90","upper.90","lower.95","upper.95","change")
data.matrix(b)
mat[,1] <- b[,1]
mat[,2] <- b[,2]
mat[,3] <- b[,3]
mat[,4] <- b[,4]
mat[,5] <- b[,5]
mat[,6] <- b[,6]

## graphic parameters
xlim <- c(-.05, .1)    # x-axis range
xlabd <- c("-.05", "0", ".05", ".1")
ylim <- c(9.5,.5)  # y-axis range
lwd=1.5
cexp=1
ptype <- 19

reg.names <- c("Surplus Government (0->1)", "Minority Government (0->1)", "Government Parties (1->4)", "Time Left in CIEP (1.0->0.5)", "Government Tenure (12->24)", "Time Since NCM (12->36)", "Number of NCMs (3->1)", "GDP P.C. Growth (2->5)","Effective Parties (2.4->4)")
change <- c("+46.3%", "+42.1%", "-29.1%","+257.5%", "-48.1%", "-36.3%", "-3.5%", "-14.9%", "-12.2%")

par(fig=c(0,1,0,1),mar=c(3.5,10.25,2,4))
plot(0,0,type="n",xlim=xlim,ylim=ylim,axes=F,ann=F,cex=cexp)

### shades
shadeColor <- c("gray85", "gray92")
xrng <- c(xlim[1],xlim[2],xlim[2],xlim[1])
polygon(y=c(0.5,0.5,1.5,1.5), x=xrng,col= shadeColor[1], border=F)
polygon(y=c(1.5,1.5,2.5,2.5), x=xrng,col= shadeColor[2], border=F)
polygon(y=c(2.5,2.5,3.5,3.5), x=xrng,col= shadeColor[1], border=F)
polygon(y=c(3.5,3.5,4.5,4.5), x=xrng,col= shadeColor[2], border=F)
polygon(y=c(4.5,4.5,5.5,5.5), x=xrng,col= shadeColor[1], border=F)
polygon(y=c(5.5,5.5,6.5,6.5), x=xrng,col= shadeColor[2], border=F)
polygon(y=c(6.5,6.5,7.5,7.5), x=xrng,col= shadeColor[1], border=F)
polygon(y=c(7.5,7.5,8.5,8.5), x=xrng,col= shadeColor[2], border=F)
polygon(y=c(8.5,8.5,9.5,9.5), x=xrng,col= shadeColor[1], border=F)

### vertical dotted line at 0
abline(v=0,lty="dotted",col="black")
idy <- seq(1,9)

### lines
for(i in 1:9){
  lines(x=c(mat[i,4],mat[i,5]), y=c(i,i),lwd=lwd)  
  points(x=mat[,1],y=idy,pch=ptype,bg="white",cex=1.5)
}

### axis
axis(2,labels=reg.names,at=idy,las=1,cex.axis=.8)
axis(4,labels=change,at=idy,las=1,cex.axis=.8)
axis(1,cex.axis=0.8, mgp=c(3,.4,0), labels=xlabd,at=xlabd)
axis(3,cex.axis=0.8, mgp=c(3,.6,0), labels=xlabd,at=xlabd)
mtext(expression(paste(Delta," in Pr(Propose)")), side=1,line=2,cex=cexp, las=1)


###############################################################################################
## Figure 3: Changes in the predicted probability of median legislator accepting a no-confidence motion over the baseline probability
###############################################################################################
a <- read.csv("prob_a.csv", sep=",", header=TRUE)

mata <- matrix(NA,5,ncol=6)         
colnames(mata) <- c("pr","lower.90","upper.90","lower.95","upper.95","change")
data.matrix(a)
mata[,1] <- a[,1]
mata[,2] <- a[,2]
mata[,3] <- a[,3]
mata[,4] <- a[,4]
mata[,5] <- a[,5]
mata[,6] <- a[,6]

## graphic parameters
xlim <- c(-.4, .8)    # x-axis range
xlabd <- c("-.4", "0", ".4", ".8")
ylim <- c(5.5,.5)  # y-axis range
lwd=1.5
cexp=1
ptype <- 19

reg.names <- c("Surplus Government (0->1)", "Minority Government (0->1)", "Government Parties (1->4)", "Time Left in CIEP (1.0->0.5)", "Returnability Index (.50->.85)")
change <- c("+863.8%", "+92.0%", "-84.6%", "+115.0%", "+408.7%")

par(fig=c(0,1,0,1),mar=c(3.5,10.25,2,3.9))
plot(0,0,type="n",xlim=xlim,ylim=ylim,axes=F,ann=F,cex=cexp)

### shades
shadeColor <- c("gray85", "gray92")
xrng <- c(xlim[1],xlim[2],xlim[2],xlim[1])
polygon(y=c(0.5,0.5,1.5,1.5), x=xrng,col= shadeColor[1], border=F)
polygon(y=c(1.5,1.5,2.5,2.5), x=xrng,col= shadeColor[2], border=F)
polygon(y=c(2.5,2.5,3.5,3.5), x=xrng,col= shadeColor[1], border=F)
polygon(y=c(3.5,3.5,4.5,4.5), x=xrng,col= shadeColor[2], border=F)
polygon(y=c(4.5,4.5,5.5,5.5), x=xrng,col= shadeColor[1], border=F)

### vertical dotted line at 0
abline(v=0,lty="dotted",col="black")
idy <- seq(1,5)

### lines
for(i in 1:5){
  lines(x=c(mata[i,4],mata[i,5]), y=c(i,i),lwd=lwd)  
  points(x=mata[,1],y=idy,pch=ptype,bg="white",cex=1.5)
}

### axis
axis(2,labels=reg.names,at=idy,las=1,cex.axis=.8)
axis(4,labels=change,at=idy,las=1,cex.axis=.8)
axis(1,cex.axis=0.8, mgp=c(3,.4,0), labels=xlabd,at=xlabd)
axis(3,cex.axis=0.8, mgp=c(3,.6,0), labels=xlabd,at=xlabd)
mtext(expression(paste(Delta," in Pr(Accept)")), side=1,line=2,cex=cexp, las=1)

###############################################################################################
## Figure 4: Probability of opposition proposing a no-confidence motion across government tenure for three different values of probability of accept.
###############################################################################################
pr_tenure = read.dta("pr_tenure.dta", convert.underscore=TRUE) 
attach(pr_tenure)

xlim <- c(0,48)    # x-axis range
ylim <- c(0,.04)  # y-axis range
lwd=1
cexp=0.8
xla <- c("0","24","48")
yla <- c("0", ".02", ".04")
ylabat <- c("", "", "")

# Unfavorable = Pr(A)=.05
par(fig=c(.1,.45,0.15,1),mar=c(1,4,2,1))
plot(vectaxis, b.best, xlim=xlim, ylim=ylim, cex=cexp, ann=F, yaxt="n", xaxt="n", type="n")

lines(vectaxis, b.worst, lty=1, lwd=1.5)
lines(vectaxis, l.worst, lty=2, lwd=1.5)
lines(vectaxis, u.worst, lty=2, lwd=1.5)
axis(1,cex.axis=cexp, mgp=c(3,.4,0), labels=xla,at=xla)
axis(2,labels=yla,at=yla,las=1,cex.axis=cexp)
mtext("Pr(A)=.05", 3, line = 1, las=1)
mtext("Pr(Propose)", 2, line = 2.5, las=1)

# Average = Pr(A)=.25
par(fig=c(.45,.725,0.15,1), new=TRUE, mar=c(1,1,2,1))
plot(vectaxis, b.okay, xlim=xlim, ylim=ylim, cex=cexp, ann=F, yaxt="n", xaxt="n", type="n")
lines(vectaxis, b.okay, lty=1, lwd=1.5)
lines(vectaxis, l.okay, lty=2, lwd=1.5)
lines(vectaxis, u.okay, lty=2, lwd=1.5)
axis(1,cex.axis=cexp, mgp=c(3,.4,0), labels=xla,at=xla)
mtext("Pr(A)=.25", 3, line = 1, las=1)
mtext("Government Tenure", 1, line=2, las=1)

# Favorable = Pr(A)=.5
par(fig=c(.725,1,0.15,1), new=TRUE, mar=c(1,1,2,1))
plot(vectaxis, b.worst, xlim=xlim, ylim=ylim, cex=cexp, ann=F, yaxt="n", xaxt="n", type="n")
lines(vectaxis, b.best, lty=1, lwd=1.5)
lines(vectaxis, l.best, lty=2, lwd=1.5)
lines(vectaxis, u.best, lty=2, lwd=1.5)
axis(1,cex.axis=cexp, mgp=c(3,.4,0), labels=xla,at=xla)
mtext("Pr(A)=.5", 3, line = 1, las=1)

# Legend
par(fig=c(0.1,1,0,.15), new=TRUE, mar=c(1,4,0,1))
plot(vectaxis,b.best,xlim=xlim,ylim=ylim,cex=cexp,type="n",axes=F,ann=F)
leg.txt <- c("Predicted Probability", "95% Confidence Interval")
lty.type <- c(1,2)
legend(0, 0.02, leg.txt, lty=lty.type, cex=1, bty="n", ncol=2)

detach(pr_tenure)

###############################################################################################
## Figure 5: Probability of opposition proposing a no-confidence motion across the electoral cycle for three different values of probability of accept
###############################################################################################
pr_ciep = read.dta("pr_ciep.dta", convert.underscore=TRUE) 
attach(pr_ciep)

xlim <- c(1,0)    # x-axis range
ylim <- c(0,.08)  # y-axis range
lwd=1
cexp=0.8
xla <- c("1",".5","0")
yla <- c("0", ".04", ".08")
ylabat <- c("", "", "")

# Unfavorable = Pr(A)=.05
par(fig=c(.1,.45,0.15,1),mar=c(1,4,2,1))
plot(vectaxis, b.best, xlim=xlim, ylim=ylim, cex=cexp, ann=F, yaxt="n", xaxt="n", type="n")

lines(vectaxis, b.worst, lty=1, lwd=1.5)
lines(vectaxis, l.worst, lty=2, lwd=1.5)
lines(vectaxis, u.worst, lty=2, lwd=1.5)
axis(1,cex.axis=cexp, mgp=c(3,.4,0), labels=xla,at=xla)
axis(2,labels=yla,at=yla,las=1,cex.axis=cexp)
mtext("Pr(A)=.05", 3, line = 1, las=1)
mtext("Pr(Propose)", 2, line = 2.5, las=1)

# Average = Pr(A)=.25
par(fig=c(.45,.725,0.15,1), new=TRUE, mar=c(1,1,2,1))
plot(vectaxis, b.best, xlim=xlim, ylim=ylim, cex=cexp, ann=F, yaxt="n", xaxt="n", type="n")

lines(vectaxis, b.okay, lty=1, lwd=1.5)
lines(vectaxis, l.okay, lty=2, lwd=1.5)
lines(vectaxis, u.okay, lty=2, lwd=1.5)
axis(1,cex.axis=cexp, mgp=c(3,.4,0), labels=xla,at=xla)
mtext("Pr(A)=.25", 3, line = 1, las=1)
mtext("Proportion of Time Left in CIEP", 1, line=2, las=1)

# Favorable = Pr(A)=.5
par(fig=c(.725,1,0.15,1), new=TRUE, mar=c(1,1,2,1))
plot(vectaxis, b.best, xlim=xlim, ylim=ylim, cex=cexp, ann=F, yaxt="n", xaxt="n", type="n")

lines(vectaxis, b.best, lty=1, lwd=1.5)
lines(vectaxis, l.best, lty=2, lwd=1.5)
lines(vectaxis, u.best, lty=2, lwd=1.5)
axis(1,cex.axis=cexp, mgp=c(3,.4,0), labels=xla,at=xla)
mtext("Pr(A)=.5", 3, line = 1, las=1)

# Legend
par(fig=c(0.1,1,0,.15), new=TRUE, mar=c(1,4,0,1))
plot(vectaxis,b.best,xlim=xlim,ylim=ylim,cex=cexp,type="n",axes=F,ann=F)
leg.txt <- c("Predicted Probability", "95% Confidence Interval")
lty.type <- c(1,2)
legend(1, 0.03, leg.txt, lty=lty.type, cex=1, bty="n", ncol=2)

detach(pr_ciep)

###############################################################################################
## Figure 6: The probability of opposition proposing a no-confidence motion for one government versus three governments illustrating the contrasting effects of government tenure and election cycle.
###############################################################################################
pr = read.dta("pr.dta", convert.underscore=TRUE) 
pr1 <- subset(pr, govt==1)
pr2 <- subset(pr, govt==2)
pr3 <- subset(pr, govt==3)

# 2X1 figure with shades
xlim <- c(1,0)    # x-axis range
ylim <- c(0,.12)  # y-axis range
lwd=2
cexp=0.8
xla <- c("1",".75",".5",".25","0")
xlabel <- c("t", "t+1", "t+2", "t+3", "t+4")
yla <- c("0", ".06", ".12")
ylabat <- c("", "", "")
shadeColor <- c("gray92", "gray85", "gray66")
xrng <- c(xlim[1],xlim[2],xlim[2],xlim[1])
yrng <- c(ylim[1],ylim[2],ylim[2],ylim[1])

# Same government across election cycle
par(fig=c(.1,1,0.57,1),mar=c(.25,4,.5,1))
plot(pr$vectaxis, pr$b, xlim=xlim, ylim=ylim, cex=cexp, ann=F, type="n",axes="F")
polygon(y=c(0,0,.12,.12), x=xrng,col= shadeColor[1], border=F)
lines(pr$vectaxis, pr$b.1, lty=1, lwd=lwd)
lines(pr$vectaxis, pr$l.1, lty=2, lwd=lwd)
lines(pr$vectaxis, pr$u.1, lty=2, lwd=lwd)
axis(2,labels=yla,at=yla,las=1,cex.axis=cexp)
mtext("Pr(Propose)", 2, line = 2.5, las=1)

# 3 governments across election cycle
par(fig=c(.1,1,0.15,.57), new=TRUE, mar=c(1,4,.25,1))
plot(pr$vectaxis, pr$b, xlim=xlim, ylim=ylim, cex=cexp, ann=F, type="n",axes="F")
polygon(x=c(1,1,.5,.5), y=yrng,col= shadeColor[1], border=F)
polygon(x=c(.5,.5,.3,.3), y=yrng,col= shadeColor[2], border=F)
polygon(x=c(.3,.3,0,0), y=yrng,col= shadeColor[3], border=F)
lines(pr1$vectaxis, pr1$b, lty=1, lwd=lwd)
lines(pr1$vectaxis, pr1$l, lty=2, lwd=lwd)
lines(pr1$vectaxis, pr1$u, lty=2, lwd=lwd)
lines(pr2$vectaxis, pr2$b, lty=1, lwd=lwd)
lines(pr2$vectaxis, pr2$l, lty=2, lwd=lwd)
lines(pr2$vectaxis, pr2$u, lty=2, lwd=lwd)
lines(pr3$vectaxis, pr3$b, lty=1, lwd=lwd)
lines(pr3$vectaxis, pr3$l, lty=2, lwd=lwd)
lines(pr3$vectaxis, pr3$u, lty=2, lwd=lwd)
axis(1,cex.axis=cexp, mgp=c(3,.4,0), labels=xlabel,at=xla)
axis(2,labels=yla,at=yla,las=1,cex.axis=cexp)
mtext("Pr(Propose)", 2, line = 2.5, las=1)

# Legend
par(fig=c(0.1,1,0,.17), new=TRUE, mar=c(1,4,1,1))
plot(pr$vectaxis, pr$b, xlim=xlim,ylim=ylim,cex=cexp,type="n",axes=F,ann=F)
leg.txt <- c("Predicted Probability", "95% Confidence Interval")
lty.type <- c(1,2)
legend(1, 0.05, leg.txt, lty=lty.type, cex=1, bty="n", ncol=2, lwd=lwd)
mtext("Election Cycle", 3, line = -1, las=1)




